import pandas as pd
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# 读取数据，假设文件名为 '附件一（训练集）.xlsx'
data = pd.read_excel('附件一（训练集）.xlsx')
print(data.columns)
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用SimHei字体
plt.rcParams['axes.unicode_minus'] = False  # 正确显示负号
# 筛选正弦波形的数据
sine_wave_data = data[data['励磁波形'] == '正弦波']

# 提取需要的变量
# 提取需要的变量
temperature = sine_wave_data['温度，oC'].values   # 温度
frequency = sine_wave_data['频率，Hz'].values     # 频率
core_loss = sine_wave_data['磁芯损耗，w/m3'].values  # 磁芯损耗
Bm = sine_wave_data.iloc[:, 4:].max(axis=1).values  # 磁通密度峰值

#比赛结束前最后一天售后群发布无水印可视化结果+无标注代码【可直接提交】为了防止倒卖， 论文写作过程中遗留数个致命问题，无关代码，该问题解决方式仅在官网授权售后群答疑，盗卖方式购买资料不提供答疑，感谢理解2024年研赛资料助攻购买链接+说明https://docs.qq.com/doc/p/e03d6bffc00b392932b7dc1fbcf746ef9b086a8a
# 拟合传统斯坦麦茨方程 P = k1 * f^a * Bm^b
def steinmetz_fun(input, k1, a, b):
    f, Bm = input
    return k1 * f**a * Bm**b

# 使用 curve_fit 进行拟合
initial_params = [1, 1.5, 2.5]  # 初始值 [k1, a, b]
input_data = np.vstack((frequency, Bm)).T
x_traditional, _ = curve_fit(steinmetz_fun, input_data.T, core_loss, p0=initial_params)

# 提取拟合的系数
k1, a, b = x_traditional

print(f'传统斯坦麦茨方程拟合参数: k1 = {k1:.4f}, a = {a:.4f}, b = {b:.4f}')

# 构造修正后的斯坦麦茨方程 P = k1 * f^a * Bm^b * (1 + c * T)
def corrected_fun(input, k1, a, b, c):
    f, Bm, T = input
    return k1 * f**a * Bm**b * (1 + c * T)

# 使用 curve_fit 进行拟合，增加温度修正项
initial_params_corrected = [k1, a, b, 0.01]  # 初始值 [k1, a, b, c]
input_data_corrected = np.vstack((frequency, Bm, temperature)).T
x_corrected, _ = curve_fit(corrected_fun, input_data_corrected.T, core_loss, p0=initial_params_corrected)

# 提取修正方程的系数
k1_corrected, a_corrected, b_corrected, c_corrected = x_corrected
#比赛结束前最后一天售后群发布无水印可视化结果+无标注代码【可直接提交】为了防止倒卖， 论文写作过程中遗留数个致命问题，无关代码，该问题解决方式仅在官网授权售后群答疑，盗卖方式购买资料不提供答疑，感谢理解2024年研赛资料助攻购买链接+说明https://docs.qq.com/doc/p/e03d6bffc00b392932b7dc1fbcf746ef9b086a8a
print(f'修正斯坦麦茨方程拟合参数: k1 = {k1_corrected:.4f}, a = {a_corrected:.4f}, b = {b_corrected:.4f}, c = {c_corrected:.4f}')

# 计算传统方程的预测值
pred_traditional = steinmetz_fun(input_data.T, *x_traditional)

# 计算修正方程的预测值
pred_corrected = corrected_fun(input_data_corrected.T, *x_corrected)

# 计算两种模型的均方误差
error_traditional = np.mean((core_loss - pred_traditional)**2)
error_corrected = np.mean((core_loss - pred_corrected)**2)

print(f'传统斯坦麦茨方程的均方误差: {error_traditional:.4f}')
print(f'修正斯坦麦茨方程的均方误差: {error_corrected:.4f}')

# 绘制实际值与预测值的对比图
plt.figure()
plt.plot(core_loss, 'o-', label='实际磁芯损耗')
plt.plot(pred_traditional, 'x-', label='传统斯坦麦茨方程预测')#比赛结束前最后一天售后群发布无水印可视化结果+无标注代码【可直接提交】为了防止倒卖， 论文写作过程中遗留数个致命问题，无关代码，该问题解决方式仅在官网授权售后群答疑，盗卖方式购买资料不提供答疑，感谢理解2024年研赛资料助攻购买链接+说明https://docs.qq.com/doc/p/e03d6bffc00b392932b7dc1fbcf746ef9b086a8a
plt.plot(pred_corrected, 's-', label='修正斯坦麦茨方程预测')
plt.legend()
plt.xlabel('样本编号')
plt.ylabel('磁芯损耗 (W/m^3)')
plt.title('实际磁芯损耗与预测磁芯损耗的对比')
plt.grid()
plt.show()
